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Abstract 

In recent years, compressed sensing (CS) has been applied in the field of synthetic aperture radar 
(SAR) imaging and shows great potential. The existing models are, however, based on application of the 
sensing matrix acquired by the exact observation functions. As a result, the corresponding reconstruction 
algorithms are much more time consuming than traditional matched filter (MF) based focusing methods, 
especially in high resolution and wide swath systems. In this paper, we formulate a new CS-SAR imaging 
model based on the use of the approximated SAR observation deducted from the inverse of focusing 
procedures. We incorporate CS and MF within an sparse regularization framework that is then solved 
by a fast iterative thresholding algorithm. The proposed model forms a new CS-SAR imaging method 
that can be applied to high-quality and high-resolution imaging under sub-Nyquist rate sampling, while 
saving the computational cost substantially both in time and memory. Simulations and real SAR data 
applications support that the proposed method can perform SAR imaging effectively and efficiently 
under Nyquist rate, especially for large scale applications. 

Index Terms 

Synthetic Aperture Radar; Compressed Sensing; Matched Filtering; Approximated Observation. 

I. Introduction 

Synthetic aperture radar (SAR) is an active microwave radar which can achieve high-resolution 
images in all time of day and weather In a SAR system, the radar emits a sequence of pulses 
along its path and receives the echoes (raw data) scattered from the targets. The reconstruction 
of the scene is traditionally achieved by matched filter (MF) based focusing algorithms which 
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are efficient but need Nyquist rate samples of the echoes. The SAR imaging with increasing 
resolution and swath requires more and more measurements, storage and downlink bandwidth. 
The current system hardware, however, frequently hampers such high-dimensional application. 

The recent development of compressed sensing (CS) brings possibility of reconstructing sparse 
or compressible signals with fewer measurements than that Nyquist requires [I3-[l4l. Several 
applications on radar system appear in recent years, which primarily concern how the data 
acquisition way can be simplified by using CS 0] [[6l] and what the potential applications will 
renovate radar imaging with CS technique [jTl flSl. Further, in the study of CS-SAR, much 
attention has been paid to the effective use of the specific SAR geography and signal form, say, 
in JH, a SAR raw data compression framework based on CS was suggested by sampling the data 
in frequency domain. An extension of this work was given in [lOJ by using the fact that very 
bright objects are always sparse, resulting in a hybrid sparse model. These works, however, do not 
apply to the CS-SAR system practically where sampling is expected in time-domain. In [fTlll . CS 
was applied on azimuth after the range compression. By combining range MF, the method was 
much more efficient, while, the redundant information in range has not been effectively utilized. 
More general CS-SAR model were reported in [1T2) by discretizing the SAR observation 
function exactly into an observation matrix, while solving by CS straightforwardly. 

All those works strongly demonstrated that some exclusive advantages of CS-SAR do exist as 
compared to the traditional SAR imaging methodologies, say, relaxation of required measure- 
ments, reduction of side lobe and a further suppression of noise [fT4ll . However, in all applications, 
a serious drawback has been observed: as compared to the traditional MF based methods, the 
computational complexity and memory cost of the CS-SAR models are much higher, so that it 
is very inefficient to be applied to high-dimensional applications. 

In this paper, we formulate a new CS-SAR framework within which the computational com- 
plexity of the CS-SAR imaging can be significantly reduced. Our main idea is to replace the 
exact observation function in the CS-SAR framework with approximated observations derived 
from the inverse of traditional MF based procedures. Such inversion has ever been applied to 
yield raw signals (the echoes) in a more economical way [fTSll [fT6ll . but requires high accuracy 
of the adopted method. In this paper, we take a further step by incorporating it into the CS 
framework, which demands only a well focusing ability to ensure CS reconstruction. We propose 
to implement the CS-SAR imaging through the sparse regularization scheme which is then solved 
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by an iterative thresholding algorithm (ITA). Accordingly, the fast speed and high efficiency of 
the new method are guaranteed respectively from the use of the approximated observation and 
from thes CS reconstruction procedure. We show that the new CS-SAR imaging method can not 
only acquire high-quality and high-resolution images with significantly reduced measurements, 
but also reduces the memory cost to 0{n) and computational complexity of one-step iteration 
to (9(nlogn), achieving the same order with the traditional SAR imaging methods. 

The reminder of the paper is organized as follows. In Section 2, we introduce the background 
knowledge on the stripmap mode SAR system and the classical CS-SAR model. In section 3, 
we present the approximated observation by calculating the inverse of MF imaging procedure. In 
Section 4, we formulate the new CS-SAR imaging method through hybridizing the approximated 
observation and sparse regularization. In Section 5, we show the simulation and application results 
of the suggested method. Conclusions are then presented in Section 6 with some useful remarks. 

Notation: We will use the subsequent notations throughout the paper: Column vectors, matrices 
and operators will be denoted respectively by bold lower case, x, bold upper case. A, and roman 
upper case, C. A^,A*,A^ denotes the transpose, conjugate and Hermitian transpose of A, 
respectively. 

II. CS-SAR MODELS BASED ON EXACT OBSERVATION 

In this section, some preliminary knowledge of CS-SAR imaging is summarized. We focus on 
the general formalization of CS-SAR model, with a more detailed introduction of the iterative 
thresholding procedures for solution of the CS-SAR models. 

A. Stripmap Mode SAR Model 

In the stripmap mode SAR, the antenna is pointed to a fixed direction and the platform flights 
with constant velocity v. Then, a complex baseband Pc{t), usually chirp, is modulated to real 
pulse p{t) = cos(27r/or + 0(r))(— ^ ^ t < ^) (/o is the carrier frequency, r is the range 
time, Wr is the elevation weight and tg is the pulse duration) and transmitted at a constant 
pulse repetition frequency (PRF). The received backscattered energy can then be modeled as a 
convolution of the pulse waveform with the ground reflectivity function, given by ifTTll 

s{r], t) = Wr{r)a{r], r) ® h{r], r) + no{r], r) (1) 
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where rj.r are respectively the azimuth and range time, no denotes the additive noise, h{r],T) is 
the time- variant convolution kernel which can be composed as: 

h{r], t) = hi{r], r) ® /i2(r/, r) (2) 

In (O, hi{r],T) is the two-dimensional azimuth modulation which is responsible for the along- 
track observation while /i2(^,t) is range convolution kernel that is identical to the transmitted 
pulses. 

Further, we can sample the continuous-time analog echo 5(77, r) and discrete the reflectivity 
map cr(7y,r), into two-dimensional arrays Y G C^^^^^ and X G C^^^^^. And then we obtain 
the following observation model for the strip mode SAR: 

y = Hx + no (3) 

where y = vec(Y) G C^^-^, / = x n!^, x = vec(X) G C^^-^, n = x n^, H is the observation 
matrix acquired from the discrete weight of ([D (more detailed information and construction of 
the observation model can be seen in [fTTl ifTSl ). and no is the noise. 

B. Formulation of CS-SAR models 

In a CS-SAR model, the data y is sampled and compressed with a proper sampling matrix 
G W^^\ m ^ n, resulting in 

y, = 0Hx + n, (4) 

When X is a sparse signal, say, most of the entries of x are zeros, the theory of CS tells when 
and how it can be recovered from the above undetermined linear system with fewer measurements 
than Nyquist criterion requires [Q 01. Generally, considering an ill-posed linear system y^ = Ax 
(A = 0H) where x is sparse enough, if the sensing matrix A satisfies some conditions like RIP 
lfT9il . X can be exactly recovered from ys with the Lq (quasi-norm) (0 < q < 1) optimization: 

min ||x||g s.t. y^ = Ax (5) 

X 

To solve ([5]), we usually use an equivalent regularization scheme with the following optimization 
problem 

min {||y,-Ax||^ + A||x||^} (6) 



February 14, 2013 



DRAFT 



5 



where A is a regularization parameter. The optimization can be efficiently solved by iterative 
thresholding algorithm (ITA) [l20l -[l22]|. In detail, an ITA generates a sequence of approximates 
according to: 

x(^+i) = E,,,^(x« + M"(y. - Ax«)) (7) 

where /x is a normalized parameter which controUs the convergence of the iteration. In d?]), 
(a = A/i) is a so-called thresholding operator which is componentwisely defined by 

E^,^(x) = (e^,^(xi), eg,a(x2), ...eg,a(x^))^ (8) 

where Cq^o- can be analytically specified when q = 0, For example, the widely used 

soft-thresholding, which corresponds to g = 1, is 

{sgn(x)(|x| — a), if \x\ > a 
(9) 
0, otherwise 

The iteration (|7]) with ([8]) is the fundamental procedure we suggest to use for the CS-SAR 
imaging. 

It can be seen that the main computation load in implementation of (|7]) comes from the 
calculation of time domain correlation A^Ax. From the viewpoint of SAR signal processing, 
this corresponds to the backscattered projection procedure, which is known to be inefficiency 
for reconstruction, even implemented by convolution as in ([I]). On the other hand, we notice 
that there exists efficient focusing methods using MF in traditional SAR signal processing. This 
type of processing is in the frequency domain, which is much faster. Moreover, unlike the 
MF based method which is usually decoupled, the sensing matrix H in (|7]), owns intrinsically 
two-dimensional structure that has to be collected and stored before imaging. Although some 
compression can be incorporated according to the structure of the matrix, it still consumes a 
huge memory. All these difficulties then hamper effective applications of the known CS-SAR 
imaging. 

The aim of the present research is to suggest a new CS-SAR imaging method, which replaces 
the use of the exact observation model by an well-defined approximation, and then makes it 
possible to reconstruct the sparse scene x via a sequence of 1-D operations. Thus, the very high 
cost of calculation and memory of the existing CS-SAR imaging methods can be significantly 
reduced. 
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III. Approximated Observation 

In this section, we first explain why an approximated observation is needed and feasible, and 
then, we provide an example to show how an approximated observation operator can be explicitly 
constructed by virtue of a concrete example from the inverse of Range-Doppler Algorithm 
(RDA). A relation between the constructed approximation and the corresponding focus method 
is analyzed, which then serves as the basis of the development of new method in the next section. 

A. Why approximation needed 

It is known that MF is fast with 0{n\ogn) complexity, which is, among the others, mainly 
due to frequency domain operations. More precisely, if we denote M the imaging procedure by 
MF, like RDA, the SAR raw data can be well focused in some conditions by 

i = My (10) 

where M is the traditional MF imaging procedure that can be calculated through decoupling it 
into a series of 1-D operators in the frequency domain. This normally leads to an 0{n\ogn) 
complexity when fast fourier transformation (FFT) type operations are employed. 

Observing these advantages, the purpose of this paper is to accelerate the known CS-SAR 
imaging procedures, so as to achieve a comparable (at the same order) complexity with the 
traditional MF based methods. 

A natural consideration is then to look for the possibility of integrating CS and MF. However, 
a direct application on the decoupling of H is impossible, because H intrinsically possesses 2-D 
structure. Nevertheless, it is known from (fTOl) that x always approximates x, say, MH ?^ I, and 
hence, M~^, whenever exists, approximates the observation H. In viewing that M is decoupled, 
it can then naturally be expected that M"^ is decoupled, too. Thus, we can expect that under 
certain conditions, some types of approximations of H can be decoupled, so as to bring an 
0{n\ogn) complexity. 

This is why we would like to approximate the observation, and in the following, we will 
introduce the details on how to construct and what constraints an appropriate approximation 
observation. 
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Fig. 1. The relations between exact observation and approximated observation. 



B. How to construct an admissible approximated observation 

Fig. [H draws the main relations between CS-SAR observation and MF reconstructions. It can 
be seen that whenever the imaging procedure M is accurate enough, M"^ can be viewed as an 
admissible substitute of H. This provides a general principle of how the observation H can be 
remodeled and approximated by any high precision imaging (or reconstruction) procedure. We 
formalize this principle further as 

G = M-^?^H (11) 

where G is any generalized right inverse of M, and M is any a high precision imaging procedure. 
We call G henceforth an approximated observation. 

However, since there are many well known imaging procedures that provided various tradeoffs 
on imaging accuracy and complexity. We need therefore further to define the extent of accuracy 
and identify the constraints under the CS-SAR framework. To see this, let us compare the exact 
observation model and the approximated observation model 

y, = 0Hx = 0Gi (12) 

It can be seen that by using approximated observation, other than reconstruct x, we actually 
reconstruct x instead, which is assumed to be an approximation of x, when it obeys to the 
following relation 

x = PeOx + Se (13) 



February 14, 2013 



DRAFT 



8 



where o denotes the Hadamard product, denote the phase error while is the error for side 
lobe, or more severely, the artifacts from unfocusing. Formally, when ([T3]) holds, there exists an 
acceptable solution with the approximated observation model. However, to find it under CS, we 
should further emphasis on a better focusing ability of M. 

As we know, a key parameter in CS-SAR, different from traditional SAR, is the sampling rate 
which measures how a SAR system benefits from CS. The least amount of samples to ensure 
the reconstruction is incoherently determined by the spars ity of scene x, usually irrelevant of 
the distribution and phases of targets in the scene. Thus, the difference of the sparsity between 
X and X determines when and how much the additional measurements does the approximated 
observation based CS-SAR methods need. It can be immediately seen from ([T3]) that the dif- 
ference is uniquely characterized via s^. That is to say, whenever the side lobes reconstructed 
from M is low enough, can be ignored, then x and x can keep the same sparsity. In this 
situation, the required least sampling rate of the approximated observation model equals to the 
original model. In turn, to prevent the approximated observation based CS-SAR method from 
demanding more samples, the side lobe should be as low as possible, or equivalently to say, a 
well focusing capacity should be a criterion to determine whether a specific focusing method 
can be used to construct the approximated observation. 

The above discussions tell a fact that the construction of the approximated observation is quite 
reflexible, which can be acquired straightforwardly based on well established algorithms with 
additional requirement on the focusing ability. In the next subsection, we present a concrete 
example using Range-Doppler-algorithm (RDA) [|23]| to show how an admissible approximated 
observation G can be simply constructed based on this principle. 

C. A concrete example 

RDA is a very popular procedure for stripmap mode SAR imaging that is simple both in 
comprehension and in implementation. The procedure (under the low squint case) consists of 
three main steps (operations): the range compression, RCMC and azimuth compression. In a 
compact form, the imaging procedure M, operated on 2-D array, can then be expressed in the 
following 

X = M(Y) = Fy{P, o C(F,[P, o (YF0]FH)} (14) 
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where X (x = vec(X)) is the reconstructed 2D SAR image, F, respectively are the DFT 
matrix and inverse DFT matrix (in practice, they are implemented by FFT) to perform, the 
subscript ry, r denotes the direction of azimuth and range where the FFT performs along, and 
are the frequency domain matched filter operations along azimuth and range, which can be 
always defined respectively by 

F.iU r) = exp[-j7i/KJ^] P.(/.; rj) = exp[-j7i/ KJ^] (15) 

In (fTSl) . /ry,/r are the frequency along Doppler and range, Ka and are the azimuth FM 
rate and the pulse FM rate. In ([T4]) . C is the RCMC interpolation operator which is essentially 
a space-variant shift, and always approximated by the truncated sinc-kemel interpolation with 

U = C(V) as 

U(/„ ^) = E ^(A. r)smc{f - (r + Ar(/„ r)) (16) 

f 

where Ar is the migration (measurement in time) to be corrected, and U, V are the signals before 
and after RCMC, respectively. 

With the so specific operations in RDA procedure, we now can derive the inverse of M quite 
simply by taking the inverse of every sub-procedures. The details are as follows: 

i) The inverse of Fourier transformations F, F^ are known as the inverse transformations, 
which are given by F^, F. It is important to keep the throwaway consistent between the pairs. 

ii) It is known that phase multiplication is a unitary transformation, so that the inverse is 
the multiplication of the conjugate phase, P*, P*, and the Hadamard multiplication can still be 
applied in order. 

in) The inverse of C is difficult to achieve directly. In fact, C is approximated from the 
accurate RCMC defined in continuous range time domain. Because the trajectories of targets 
with different range gates are disjoint, this shift is a one-to-one mapping and the inverse of the 
origin RCMC exists. We can also approximate through interpolation V = D(U) that 

r) = J2 ^(U r)smcif - (r + Ar(/„ r))) (17) 

r 

Based on the above exposition, the approximated observation G deduced from RDA can then 
be explicitly expressed by: 

G(X) = {P; o (FjDp; o (F,X)]F.)}Fj^ (18) 
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We show that the so constructed approximated observation G has an interesting property: It 
is still a linear operator and its conjugate transposition equals to 

Theorem 1. G is a linear operator with the property G^ = M. 

Proof: The linearity of G and M is obvious because all the sub-operations are linear. Let 
X denote the vector form of X, namely, x = vec(X). Then, by definition, the linear operators 
G and M can be written as matrices, and we then have 

vec(G(X)) = Gx = Ff P;Ff F,DP;F^x (19) 

vec(M(Y)) =My = Ff P^CF^Ff P^F^x (20) 

where 

F^ = I,^®F^,F, = F^®I,^ (21) 

P^ = diag(vec(P^)), P, = diag(vec(P,)) (22) 

C and D are real matrices defined by 

C(^,J•) = smc(^ + Ar (.)//,) 

^/ j) mod = (23) 

D(^,i) = sinc(^ + ArO•)//.) 

In ([23]) . Ar = vec(AR), AR is the discretion of Ar. Observing from (l23l) and fs is the pulse 
sampling interval, it is easy to check that D = C^. Consequently, comparing (1201) and ([T9]) . we 
conclude that G = M^. ■ 
Theorem 1 shows that we have actually taken the conjugate transposition of M as an approx- 
imated observation of H. Such coincidence plays an important role in the new method to be 
suggested in the next section. 



D. Generalization 

The approach we have applied to define the approximated observation by the inverse of RDA 
procedure can be generalized in the following two folds: 

1) For high squint cases, we can incorporate secondary compression in RDA or derive 
an approximated observation from the inverse of other focusing methods like Chirp-Scaling 

^It is to say G is nearly unitary, since only a minor approximation on calculation of the inverse is included 
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Algorithm (CSA), uj — k algorithms, so as to enhance the focusing ability. With similar sub- 
operations, like FFT, phase multiplication and interpolation, the acquisition of the inverse is just 
as same as RDA. 

This principle can be applied to yield more general algorithms, however, we will not enumerate 
all possible extensions but remind that the decoupled structure of MF makes the inverse always 
achievable. This is the reason why MF is fast and why we propose to apply the approximated 
observation instead of the exact observation in the CS-SAR imaging system. 

2) In the above derivation, we have assumed that the transmitted signal form is standard chirp 
and M focus both the azimuth and range direction. In fact, the azimuth modulation is the exclu- 
sive property and main difficult of SAR signal processing while the range convolution, which 
possesses a simple 1-D structure, can be modeled directly. So, we can apply the approximated 
observation to non-chirp cases, by replacing P* with the transmitted pulse, which is always 
recorded in modem SAR systems. This extension is also of specifical necessity in CS-SAR 
because the design of the pulse form is also a very important issue. 

IV. CS-SAR IMAGING BASED ON APPROXIMATED OBSERVATION 

In this section, we formulate the new CS-SAR imaging method based on the use of ap- 
proximated observation. An Lq regularization model together with the fast iterative thresholding 
algorithm will be suggested. 

A. The new CS-SAR method 

By replacing the exact observation H using the approximated observation G in (fTSl) . we can 
acquire the following CS-SAR model: 

min{||Y, - 0^G(X)0,||| + A||X||^} (24) 

where || • is the Frobenius norm of a matrix, 0^ and 0^ are the sampling operators in azimuth 
and range directions, which corresponds to the general sampling operator in (|4]). It is well 
defined because the azimuth signal is of discrete form and the range signal is of continuous form, 
and thus the sampling procedure of the two types of signals are usually physically separated^. 

^Although the sampling scheme on range may vary pulse-by-pulse, we still use this expression which is easily understood 
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Then, due to the hnearity of G, the model (|24|) can still be very fast solved by ITA, which 
reads in this case that 



X(^+^) = Ei,,^(X« + /xM(0j(Y, - 0,G(X»)e.)0^)) 



(25) 



In this paper, we simply select q = 1 while parameters /x, A will be preset according to the next 
subsection. 

Fig. [21 below shows the flow chart of algorithm (l25l) . which tells that ITA provides an intuitive 
explanation in terms of SAR signal processing. It is seen that at each iteration, the ITA can be 
decomposed into mainly three procedures: the compressed data simulation, the matched filter on 
the residual and the thresholding for new estimation. Physically, this means that in every iteration 
of the ITA, the useful information in the residue (not the raw data) is first extracted by MF and 
then added to the current estimate to yield a new update, finally the thresholding procedure 
enforces the sparsity through regularizing the noise and ambiguity from under-sampling. 
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Fig. 2. An explanation of the proposed algorithm in one-step iteration. The compressed data is alternatively processed by MF 
and thresholding. 



The algorithm stops when converges or achieves the maximum number of iterations. For 
convenience of use, we list pseudo-code of ITA ([25]) as Algorithm 1 below. We further show 
how the parameters in (l25l) can be set adaptively. 

B. Parameter Setting 

There are two parameters /x and A in ([251) that need to be set. First, ji controls convergence 
of the ITA that the inverse should obeys. 

< /i"^ < ||A||^ (26) 
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Algorithm 1 : Iterative thresholding algorithm for approximate observation based CS-SAR 
imaging 

Require: SAR raw echoes Y^, approximated observation operator G and M, sampling operator 0^, 0^- 

Ensure: The recovery image X* 

Initial: X*^^^ = 0,A, /i and max iteration /max 

1: for z = to /max do 

2: Residue: R^'^ = Y, - 0^G(X^^^)0^ 

3: MF on residue: AX^^^ = M(0^R^^)0^) 

4: Thresholding: X^^+^) = Ei,a;.(X(^) + fiAX^'^) 

5: end for 



However, it is difficult to calculate || 0^00^11 2 directly, where operator G is included. As an 
alternative, we adopt the adaptive step selection strategy in [l24ll as: 

/i, = ||0,G(AX«)0.|||/||AX«|||. (27) 

where AX^^^ equals to AX^^^ at the support of X^^"^^ and equal to zero elsewhere. It is easy 
to demonstrate that /x^ satisfies (IZTl) . and as reported in [l24]| . such choice has an additional 
advantage of accelerating the algorithm. 

Further, the regularization parameter A, which functions to compromise the reconstruction 
precision and the sparsity of the solutions obtained, has a substantial impact on the imaging 
result. Fortunately, as a part of the Lq regularization theory, the optimally A has been resolved in 
|l2T]l . whenever the problem's sparsity is known. More precisely, assume the considered problem 
has sparsity k (i.e. a fc-sparsity problem), then the optimal setting problem of parameter A* is 
shown to satisfy 

A*e [|6^(x*)|,^^^,|6^(x*|,)^] (28) 

where 6/^(x) = x + /xAx, |6^(x*)|^ is its k-th largest component in magnitude. 

Therefore, we suggest the setting that in the nth iteration A^ = |6^(x(^))|^^^ //x^ ( A^/x^ is 
independent of /x^). The sparsity fc, which determines A^, can be much flexible to be set, say, 
based on a prior upper estimation on sparsity of the target scene. 
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C. Computation Cost 

Let us compare the computational complexity and the memory occupation of the suggested 
CS-SAR imaging method ([24h . as compared with the known CS-SAR model The purpose 
is to see how much reduction of computational cost of the new model has been brought. In the 
calculation, we have used some standard notations which are: the number of required iteration /, 
the sampling rate s, the number of range gates n^, the number of range lines (n = rir x Uj^), 
the number of samples of sent pulse u^, the number of samples of the synthetic aperture time 
(they equal to the time bandwidth product (TBP) in each direction), and the TBP of radar 
signal u = X Ur. 

With these notations, we can calculate the computational complexity of the approximate 
observation based CS-SAR, C^, and the computational complexity of the exact observation 
based CS-SAR, Ce as follows. For C^, it includes calculations of an inverse MF procedure and 
a MF procedure, which has commonly the computational complexity of (9(nlog2n), together 
with a decoupled thresholding operator with complexity 0{n) in a single step. Thus, the total 
cost is at the order Ca = 0(/nlog2 n). For Cg, it includes calculations of a single iteration, two 
matrix multiplications and the thresholding procedure. Since there are only few non-zero entries 
in H, say, nearly us in every column, when coding it using two-dimensional convolution, it 
needs at least 2uns complex multiplication. Thus, we find that the total cost is Ce = 0{luns). 
Then, the ratio between Ce and Co is given by re'- 

rc = 0{us/log^n) (29) 

It is seen from (1291) that the ratio rc depends linearly on the TBP of radar signal u. In SAR 
applications, the u is always designed very large (thousands even millions) to improve the 
reconstruction signal to noise ratio (SNR), which will bring very high computational cost of the 
time domain reconstruction method. 

The memory loads of the approximate observation based CS-SAR, M^, and the memory loads 
of the exact observation based CS-SAR, Me can be estimated in the subsequent way. For M^, it 
contains only the storage of input, output and several parameter matrices (i.e., azimuth matched 
filter, range matched filter and the amount of migration in RDA), which is summed up to 0{n) 
bytes memory occupation. For Mg, although no filters are stored, it needs additionally to store 
a sensing matrix, with the number of non-zero entries of uns. But, because the Doppler history 
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with same range cell share the identical patterns, we only need to store an intact holistic pattern 
(the convolution kernel) for each range gate to achieve a compression, resulting in an additional 
memory occupation of Wurir bytes (a complex number occupies 16 bytes memory), as compared 
with Ma. This additional cost can be very large in spacebome SAR systems. For example, when 
u = 10^ and = 10^, it requires more than 100 GB memory to store the array. But the memory 
cost of RDA is only a few hundreds MB in the same condition. This will further hamper the 
application of time domain methods into practice. 

Finally, the required number of iteration steps is difficult to compare analytically, but in 
practice, no obvious difference is observed. 

D. A Summary 

From the analysis in the previous subsections, we can see that the suggested new model 
(l24l) and method (125)) have constituted a more efficient CS based SAR imaging method. While 
preserving CS features, the new method has the following exclusive advantages: 

• Lower computational cost: Due to the use of approximated observation, the method only 
involves 1-D operations which makes the imaging process extremely efficient. It has reduced the 
computational complexity of the existing exact observation based method significantly, as shown 
in (1291) . Meanwhile, taking full advantages of the decoupled structure after approximating the 
observation, the proposed method can save the memory cost with a remarkable amount, which 
is sometimes of more significances. 

• Higher Compatibility: As compared with the traditional MF based SAR imaging procedure, 
the new method uses the same or similar operations, except an additional thresholding operation 
that yields the sparsity of solution. In particular, the new procedure can be seen as a successive 
iterative refinement of the well known MF based method, which makes the new method more 
consistent. As a result, the proposed model requires little modification of the existing SAR 
imaging algorithms, which makes the combination of MF and CS much simpler. 

All these features make the suggested new CS-SAR imaging method more useful and efficient, 
and particularly possible to be applied in high dimensional SAR applications. We will provide 
simulations and applications in the next section to further support such benefit. 
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V. Simulations and applications 

In this section, several simulations and applications are provided to demonstrate the effective- 
ness and efficiency of the proposed CS-SAR imaging method. For abbreviations, we denote by 
CSRDA the CS-SAR imaging method ([24l) with the approximated observation G acquired from 
the inverse of RDA, and by CSEO the CS-SAR method © with the exact observation H. 

We first conduct a series of simulations to compare the performance of the CSRDA method, the 
CSEO method and the traditional RDA method in terms of reconstruction ability, reconstruction 
quality and reconstruction cost. Then, we apply the CSRDA to some real SAR imaging tasks from 
RADARSAT-1, which then further demonstrates the outperformance of the suggested method. 

The sampling scheme used in the simulations are specified as follows. In the azimuth direction, 
we employed random downsampling, realized by selecting random rows from the raw data Y, 
with sampling rate 5^. In the range direction, we picked up random samples independently on 
each sampled echoes in azimuth, with sampling rate s^. And we keep the ratio between Sa and 
5^ as 1 : 5^. 

Table 2 lists the primary SAR parameters used in both simulations and applications. All the 
experiments were conducted on a work station of 8-core 2.4GHz CPU with 32G memory. The 
CSRDA was implemented in MATLAB 2012a while the CSEO using optimized convolution was 
implemented in C++ with parallel codes and careful array operations. 

A. Simulations 

In the simulations, the scene was taken as 180 x 180 while the scattered coefficients were 
chosen with unit amplitude and uniform random phases. The raw data were first generated in 
time-domain by exact slant range and then sampled with different rate, to yield the compressed 
measurements. The sparsity parameter k was kept the same for CSRDA and CSEO and the 
maximum iteration steps was set to 100 for both methods. The aims of simulations is then 
to compare the reconstruction ability (RA), reconstruction quality (RQ) and reconstruction cost 
(RC), of each competitive SAR imaging methods. These are measured respectively by the lowest 

^The suggested sampling strategy was designed to comprehensively compare the reconstruction algorithms. The proposed 
model itself is adaptively to more complicated sampling schemes, for example, jitter sampling in the azimuth direction I.12J and 
random demodulation 1251 in the range direction. 
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TABLE I 

Primary Parameter of S AR system and geometry 



Parameter 


Symbol 


Simulation 


RadarSat-1 


Slant range of scene center(km) 


Rc 


20 


1016.7 


Eiiective radar velocity(m/s) 


V 


350 


7062 


Beam squint angle(rad) 


a 
U 


U 


U.Ub 


Radar center frequency (MHz) 


fo 


5000 


5300 


Pulse repetition frequency(Hz) 


Fa 


175 


1256.98 


Range FM rate(MHz//is) 


Kr 


37.5 


0.72135 


Pulse duration(/is) 


Tr 


2 


41.75 


Sampling rate(MHz) 


Fr 


75 


32.317 



amount of measurements a method can successfully reconstruct an image, the side lobe and 
resolution of reconstructed point target, and the computation time cost by a method to recover 
the image. 

1) RA comparison: A set of simulations was made where 9 targets were located at the center 
with intervals of 6 samples. We varied the sampling rate ranged from 100% to 0.6% and added 
gaussian noise with level of 20 dB. We applied RDA, CSRDA and CSEO to this experiment 
with sparsity parameter k = 18. Some of the simulation results are shown in Fig. |3]and |4l 

It is seen from the top row of Fig. [3] that with full samples (namely, with 100% sampling rate), 
all the methods RDA, CSEO and CSRDA can successfully recover the scene, say, the amplitude 
of the target is maintained and no false target is observed. But the reconstruction of RDA is 
with serious side lobes, which is not observed in CSRDA and CSEO. This shows the exclusive 
advantage of the sparse regularization based CS-SAR imaging methodologies, as reported in 
|fT4ll . When we reduce the sampling rate, say, 10% samples, as seen in the bottom of Fig. [3] 
RDA fails to recover the scene, while CSEO and CSRDA both can not only perfectly reconstruct 
the scene but with significantly reduced side lobes. In this case, no visible difference can be 
observed for CSEO and CSRDA. Nevertheless, when the sampling rate continues reducing as 
in Fig. 131 we found that both CSRDA and CSEO can reconstruct the image with only 0.65% of 
the samples. However, CSRDA fails with 0.55% samples while CSEO fails until the sampling 
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Fig. 3. The reconstruction results of 9 point targets simulations with different sampling rate. From left to right are the 
reconstruction results of RDA, CSRDA and CSEO, respectively. And the top row is with full samples while the bottom with 
10% samples 



rate takes 0.45%. 

All the results consistently show that benefited from sparse regularization, the approximate 
observation based CS-SAR method can reconstruct sparse scenes with far less samples than 
Nyquist rate requires. However, caused by approximation, it requires a little more samples to 
reconstruct the scene. 

2) RQ comparison: Sparse regularization was demonstrated in SAR and CS-SAR imaging 
the ability of reducing the side lobe and simultaneously improving the resolving ability [fT4| . 
Hence, we are interested in whether the enhancement is kept when approximated observation is 
included, especially when the effect of the accuracy of the observation is excluded. To illustrate 
it, we compare the reconstruction quality of RDA and CSRDA in terms of side lobe and spatial 
resolution, when successful reconstruction is achieved. The side lobe is evaluated via the peak 
side lobe ratio (PSLR), defined as the ratio of the peak intensity of the most prominent side lobe 
to the peak intensity of the main lobe, say, the smaller the PSLR, the better an algorithm. The 
spatial resolution is measured by the impulse response width (IRW), defined as the width of the 
main lobe of the impulse response, measured by 3 dB below the peak value, or the minimum 
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Fig. 4. The detailed comparison on needed least samples to reconstruct the image. The top row is the results from CSRDA, 
and the bottom row is from CSEO. From left to right are the results corresponding to 0.65%, 0.55%, 0.45% sampling rate, 
respectively. 



distance an algorithm can separate two targets, which should be also as smaller as better. We 
have perform a one point simulation with the upsampled factor 16 while the target is analyzed 
by a 16 X 16 chip centered on the peak, to yield a more detailed analysis on both the main lobe 
and the side lobe. The sampling rate is fixed (10% for CSRDA while 100% for RDA) to ensure 
a successful reconstruction, and the sparsity parameter k is set 600. Some of the simulation 
results are given in Fig. |5l 

Fig. [5] show the contours of the reconstruction results, with contour lines of -3 dB (the boundary 
of main lobe) and -13 dB (the PSLR of traditional MF output). The comparison intuitively 
shows that both the area of the main lobe and the PSLR in the side lobe reconstructed from 
CSRDA is much smaller than those reconstructed from RDA. Further in the table, details on the 
reconstruction quality are presented in azimuth and range directions. It is seen that the width of 
main lobe reconstructed from RDA is 15 samples both in range and azimuth, but for CSRDA, 
the width is only 8 samples in the two directions. For the side lobe, the PSLR reconstructed 
from RDA is -13.32 dB in range and azimuth direction, which is very obvious. But the CSRDA 
reduces the PSLR to -21.3 dB and -22.7 dB in azimuth and range, respectively. 



February 14, 2013 



DRAFT 



20 




(a) (b) 



Method 


PSLR(R) 


PSLR(A) 


IRW(R) 


IRW(A) 


RDA 


-13.32 dB 


-13.32 dB 


15 samples 


15 samples 


CSRDA 


-22.71 dB 


-21.32 dB 


8 samples 


8 samples 



Fig. 5. Contours of magnitude. The red, yellow, light blue, blue contour lines are corresponding to the value of -3 dB,-13 
dB,-23 dB and -33 dB. (a) Contours from RDA. (b) Contours from CSRDA. (c) Contours from CSEO. 

We further demonstrate the enhancement of resolution by CSRDA in another way. We add 
two point targets in the above simulation, with respectively azimuth and range interval of 12 
samples from the center. We then test whether the three targets can be separated from CSRDA. 
It is seen from the simulation result (Fig. |6ta)) that the reconstructed result of traditional RDA 
exhibits an overlapping of the main lobe, thus the targets are inseparable. By using CSRDA in 
Fig. Ob), one can clearly distinguish the locations of the three points, demonstrating an obvious 
enhancement of the resolution. 

All of the simulations in this subsection demonstrate that by combining MF with sparse 
regularization, we can reduce the side lobe and improve the resolution simultaneously to a great 
extent. Note that, these two goals have been regarded as trade-offs traditionally, if only MF is 
employed. 

3) RC comparison: Finally, we compare the CPU time takes by CSRDA and CSEO. According 
to the analysis in section IV, the computational complexity of CSRDA and CSEO depends on 
the scene size n and TBP of radar signal u. So we generated 10 examples for each set of fixed 
scene size n and TBP of radar signal u in simulation, with a constant sampling rate s 10%. 
The average computational time in a single iteration of the two methods was then recorded. The 
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Range Samples 



Azimuth Samples 



Range Samples 



(a) (b) 
Fig. 6. Superresolution ability of the proposed method, (a) Result of RDA. (b) Result of CSRDA. 



comparison results are shown in Fig. H 





Dimension of problem 



TBP of radar signal 



(a) (b) 
Fig. 7. The CPU time of CSRDA and CSEO, where (a) u is fixed as 10^. (b) p is fixed as 10^. 

As we can observe from Fig. U\ when u is fixed as 10^, the CSRDA scales very well to very 
high dimensional problems, since even with p = 10^, it only takes several second to finish an 
iteration. The CSEO is also insensitive to dimension when u is fixed. However, the CPU time of 
CSEO is consistently higher than that of CSRDA, with a ratio around 100. On the other hand, 
it is seen from Fig. |7ta) that when n is fixed to 10^, the CPU time of CSRDA is constantly 3s, 
but, the CPU time of CSEO depends linear on u, which becomes more and more costly as u 
increases. 

This RC comparison shows that the approximate observation based CSRDA is much faster 
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than the time domain method CSEO, benefited from the (9(nlog2 n) computational cost of MR 
And the computational complexity of CSEO depends linearly on both the TBP of radar pulse 
and the scene size. So when the u is large, i.e. 10^ which is commonly in spacebome cases, 
CSRDA is expected to accelerate the CS-SAR reconstruction more than thousands of times. This 
acceleration of computational time together with the shown memory saving capacity demonstrate 
the superiority of the proposed method. 

All the above simulations support that the suggested approximated observation based CS- 
SAR method, CSRDA, is capable of high quality imaging under Nyquist rates and with a 
comparable complexity as the traditional MF based imaging method. Especially, as compared 
with the exact observation based CS-SAR imaging, CSRDA preserves the features of imaging 
under Nyquist rate and reconstruction with feature enhancement, while reducing the imaging 
complexity dramatically. Such significant complexity reduction property makes the new method 
applicable to large scale imagery application (as will be demonstrated in the next subsection). 
This is, however, at the sacrifice of reconstruction quality or, equivalently, must be compensated 
with additional measurements. Thus, the new method provides a satisfying trade-off between the 
reconstruction complexity and quality. 

B. Application 

We have applied the new method, CSRDA, along with RDA and CSEO, to some real SAR 
imaging tasks. RADARSAT-1 is a famous satellite SAR launched at 1995, and the data used 
in this application were collected on June 16, 2002 with Fine Beam 2 about Vancouver region. 
The related key parameters of SAR system are as in Table. 1. 

We first applied the 3 methods to reconstruction of the region of English Bay, in which 
6 vessels are sparsely distributed, a very typical sparse scene. The scene was digitalized as 
1024 X 512 image with azimuth resolution 9 m and range resolution 6 m. We then reconstruct 
the image by the 3 imaging methods with sparsity K = 10000, and with varied sampling rates 
from 100% to 5%. Some typical results of reconstructions are shown in Fig. [8l As expected, 
the application shows a completely similar performance as that in the simulations. For example. 
Fig. [St a) exhibits that RDA can only reconstruct the image when samples are fully adopted, but 
strong side lobe is observed. When sampling rate is 20%, RDA fails with obvious ambiguities, 
but CSRDA and CSEO both can perfectly recover the image, with much reduced side lobe. Fig. 
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[8][c)(f) then show that when samphng rate is reduced to 5%, the reconstruction of CSEO is with 
shghtly higher precision than that of CSRDA, though both can still recover the targets. On the 
other hand, as listed in Fig. Ul CSRDA exhibits its dominant advantage in computational cost as 
compared with CSEO. For example, the computation time of reconstruction, when 20% samples 
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(a) (b) 
Fig. 9. Application results on RADARSAT-1. (a) RDA with full samples, (b) CSRDA with full samples. 



(a) (b) 
Fig. 10. Detailed comparison on the selected area with enlarged scale, (a) RDA. (b) CSRDA. 

are used, by CSRDA is 1 minutes, while by CSEO is about 9 hours. Also, CSEO needs to store 
the sensing matrix which occupies about 16Gb memory, while only 100Mb for CSRDA. We 
further applied the CSRDA to the large scale imaging problem together with RDA. However, the 
CSEO is not compared in this experiment because the memory cost is beyond the computational 
ability of our computer. The scene has a size 2048 x 2500 samples, which is of large scale 
but not so sparse. CSRDA can be applied in principle because sparse regularization, which is 
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adopted in CS-SAR, can be used as a feature enhanced imaging method (with suppressed side 
lobe and improved resolution), as demonstrated in the simulations. The reconstruction results 
by RDA and CSRDA (with 100% sampling rate) are shown in Fig. [9l It is seen from Fig. [9] 
that CSRDA has resulted the reconstruction with improved resolution and reduced side lobes, 
as demonstrated in the zoomed Fig. [TOl 

The applications above support that the suggested approximated observation based CS-SAR 
imaging is effective and efficient, especially applicable to high dimensional SAR imaging appli- 
cations. Such benefit clearly improves on the currently used exact observation based CS-SAR 
imaging methods. 

VI. Conclusion 

Compressed sensing (CS) has been applied to yield novel SAR imaging methodologies under 
Nyquist sampling in recent years. The resultant CS-SAR models are time domain based and using 
the exact observation, which then makes it of very high computational cost, and is difficult to 
be applied in high dimensional applications. In this paper, we have proposed an approximated 
observation and frequency domain based CS-SAR imaging method, with which the computational 
complexity can be dramatically reduced. 

The main contributions of the present work are as follows. 

i) Instead of the exact observation matrix, an operator, called the approximated observation 
is constructed to generate SAR raw data by means of inverse of any traditional MF based 
imaging procedure (like RDA). Such generic construction makes the SAR imaging approximately 
decoupled, very fast in processing, while naturally connecting the existing MF based SAR 
imaging algorithms. 

ii) Incorporating the approximated observation into CS-SAR framework, an efficient sparse 
regularization based CS-SAR model is formulated. The new model combines naturally CS and 
MF and is compatible with most existing SAR imaging methods, which needs a little modification 
of the current SAR imaging technologies. 

iii) With the use of approximated observation, an iterative thresholding algorithm is suggested 
for fast solution of the new CS-SAR model, which forms a low complexity, CS featured, new 
SAR imaging method. 
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We have tested and applied the new suggested CS-SAR method with a series of simulations 
and applications. The experiments consistently support that the new method is capable of recon- 
structing sparse scenes with far fewer measurements than Nyquist requires, yielding always a 
feature enhanced high quality imaging and bringing a speed up of reconstruction hundreds times 
as compared with the exact observation based CS-SAR methods. Due to the fast and feature 
enhancement features, the new CS-SAR method can be accepted as an efficient CS type SAR 
imaging technique, especially for high dimensional imaging applications. 

It is worthwhile, however, to remark that although significant complexity reduction (speed- 
up of reconstruction) can be brought, the use of the approximated observation requires more 
samples to reconstruct a scene. Thus, how and what extent does the approximation affect the 
reconstruction deserves a further study. Moreover, since there are many possibilities of concrete 
realizations of the approximated observation, the criterion on how to select an appropriate one 
deserved study. All those problems are under our current research. 
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